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Abstract 

A barrier penetration model has been proposed to explain the temperature 
driven spontaneous melting of the DNA oligomers into two separate single 
strands whereas the fraction of partially melted intermediate states can be 
understood on the basis of bound state solution of the effective potential 
obtained using Peyrard-Bishop formalism. The excellent agreement of the 
predictions with the recent experimental measurements (Euro. Phys. Lett. 
62, 452, 2003, Phy Rev. Lett. 91, 148101, 2003) provides strong justification 
for the proposed model. A simple relation has been obtained that correlates 
the strand stiffness with the temperature at which the spontaneous melting 
probability becomes half. This is an important outcome as the stiffness pa- 
rameter of the double stranded conformation can be directly estimated by 
knowing the mid point temperature which can be measured experimentally. 
Further, it is shown that the derivative of this probability with respect to 
temperature is quite sensitive to the parameters of the entropic barrier and 
can be used as the experimental probes to study the stacking interaction in 
detail. 

PACS numbers:87.10.+e, 03.40.Kf, 05.90.+m 
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I. INTRODUCTION 



The DNA molecule is a double-stranded biopolymer with two complementary sugar- 
phosphate chains (backbones) twisted around each other forming a right handed helix [1,2]. 
Each chain is a linear polynucleotide consisting of four kinds of bases: Adenine, Guanine 
(purines) and Cytosine, Thymine (pyrimidines) . The two chains are joined together by 
hydrogen bonds between pairs of nucleotides, A-T with two hydrogen bonds and G-C with 
three hydrogen bonds. At sufficiently high temperature, the double stranded DNA (dsDNA) 
hehx melts and the molecules separate into two single strands (ssDNA) . The experimental 
procedure consists in preparing a set of dsDNA with the same length and sequences, which 
are immersed in an aqueous solution with given physiological conditions. Then the fraction 
of broken base pairs as a function of temperature, referred to as the melting curve, is 
measured by the UV light absorption, typically at about 260 nm. Long DNA molecules 
give rise to steps in the melting curves corresponding to different regions melting at different 
temperatures whereas synthetic oligonucleotides do not show any multistep behavior (see [3] 
for a review). This denaturation process, often interpreted as a reminiscent of a discontinuous 
first order phase transition, has attracted a great deal of interest to study the physics 
of highly nonlinear biomolecular system and also it has led to several micromanipulation 
experiments [4-27]. Prom biophysics point of view, the strand separation is a key aspect of 
DNA transcription and replication [28]. During transcription, a transient bubble of ssDNA is 
formed to allow the enzyme that makes a RNA copy of the DNA sequence to access the DNA 
bases. During replication that occurs during cell division for example, the dsDNA separates 
into two, each strand then serving as a template for the synthesis of a new dsDNA. On the 
otherhand, accurate prediction of DNAs thermal denaturation is an important parameter 
for many biological techniques like polymerases chain reaction (PGR) where dsDNA are 
routinely converted to separated ssDNAs by melting at elevated temperature. Therefore, 
the study of DNA denaturation which is an important biological phenomena has become a 
subject of interdisciplinary interest. 
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The early theoretical models like Poland- Scheraga (PS) type [4,29,30] considers the DNA 
molecule as being composed of an alternating sequence of bound and denaturated (loop or 
bubble) states. Typically, a bound state is energetically favored over an unbound one while 
a denaturated segment is entropically favored over a bound one. The DNA denaturation 
can be viewed as a competition between the entropy of the denaturated loops and the bound 
state energy of the sequence. In PS type models, the segments that compose the chain are 
assumed to be noninteracting with each other. This assumption considerably simplifies the 
calculation of free energy. In the past, the entropy of the denaturated loops has been eval- 
uated by modeling them either as ideal random walks [29] or as self avoiding walks [30] . It 
has been found that within this approach the denaturation transition of DNA is continuous 
both in two and three dimensions, and it becomes first order only for d> 4. A denaturated 
segment in between two bounded structures may introduce flexibility , allowing neighboring 
segments to interact with the bubble and with each other. Inclusion of this excluded volume 
interaction effect in PS model drives the transition first order even ior d — 2 { and above) 
[15,16,23,27]. In a different approach, Peyrard and Bishop (PB) [5] proposed a microscopic 
model to describe the dynamics of the DNA denaturation by introducing an interaction 
potential which depends on the transverse stretching of the hydrogen bonds between the 
complementary base pairs. An improved version of this PB model with more realistic treat- 
ment of the backbone stiffness [6,7] yields what appears to be a first order phase transition. 
Despite the simple formalism, the PB model has successfully reproduced the essential fea- 
tures of thermally induced denaturation of long DNA chains. It has been used to estimate 
the melting curves of very short heterogeneous DNA segments, in excellent quantitative 
agreement with experimental data [11]. It also provides the characteristic multistep melting 
behavior observed in heterogeneous DNA sequence [9]. Recently, this model has been used 
to investigate charge transport properties in DNA chain [31] and the role of nonlinearity and 
entropy provided by the local base pairs on the interaction of the DNA transcription [32]. 

As mentioned before, the UV absorption of the spectroscopic measurements are inter- 
preted as the yield for the average fraction of the open base pairs. The denaturation process 
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leads to a sharp transition from dsDNA to ssDNA when the sequence is long. For short 
oligomers, this transition is not sharp enough. Due to finite size effect, two fractions, fully 
dissociated molecules (p) due to opening of all the base pairs spontaneously and partially 
melted intermediate states (a), may co-exist over a wide range of temperature below the 
melting point. Prom the spectroscopic measurements, it is not possible to estimate these 
two fractions separately. Recently, a new method has been proposed [33] to quantify the 
presence of intermediate states and fully open molecules by combining the UV spectroscopy 
with the quenching technique. The sequences which are partially self-complementary are 
prepared in the duplex form by hybridizing with the exact reverse complement. On heat- 
ing, the duplex either melts spontaneously into two separate single strands or remains as 
connected but with many broken base pairs. Molecules which are completely open form 
hairpins on quenching while the molecules which are partially open close again as duplexes. 
The relative amount of hairpins give the fraction p of completely open molecules whereas the 
UV spectroscopic measurements yield the fraction of total open base pairs /. The difference 
cr — f — p represents the fraction of partially melted intermediate states. This method has 
yielded new data on the fraction of /, p and a and also the relative average bubble size 
< Z > as a function of temperatures in the melting transition of several DNA oligomers 
[33,34]. In this work, we try to understand this new and exciting experimental data both 
qualitatively and quantitatively within the frame work of PB formalism which has been quite 
successful elsewhere. We avoid the rigorous solution of the transfer integral (TI) equation 
to evaluate the configurational partition function. Instead, we concentrate on yet simple 
but well approximated TI integral which yields a Schrodinger like second order differential 
equation. Por simplicity, we use Morse potential [35] with inclusion of anharmonic stacking 
interaction [6]. We propose here that for the spontaneous melting of the dsDNA leading to 
a complete strand separation , the molecules are required to overcome an effective barrier 
through thermally assisted tunneling whereas the fraction of partially melted intermediate 
states can be explained on the basis of the bound state solution of the effective potential 
in the temperature range for which the bound state exists. The excellent agreement of this 
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prediction with experiments [33,34] provides strong justification for the proposed formahsm. 
Since the spontaneous melting probabihty is considered as the tunnehng phenomena, there 
exists a mid point temperature at which the thermal energy becomes equal to the height of 
the entropic barrier. This aspect has been utilized further to derive a simple but an impor- 
tant relation between the double strand stiffness parameter and the mid point temperature 
which can be determined experimentally. It is shown that the height and the width of the 
derivative of the double strand separation probability with respect to temperature depend 
strongly on the parameters of the entropic barrier and thus, can be used as the experimental 
tool to learn more about the nature of the stacking interaction. 



II. FORMALISM 

The details of the PB model with various improvements have been described elsewhere 
[5-8,12,22] and also we have mentioned it briefly in the appendix for completeness. What 
concerns here in the present study is the partition function of the configurational part given 

by 

1 roo 

Zy = — nf^,ciy„ e-/^^("-), (1) 

OL J —oo 



where the constant a having the dimension of length has been introduced to make Zy 
dimensionless. The potential energy component H{yn, yn-i) is 



n 



+Dn [exp {-anVn) - 1]^ - Dn. 

(2) 



The index n runs over all the base pairs and -\/2?/n is the relative stretching of the 
hydrogen bond from the equilibrium position of the n*^ base pair. The first term represents 
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the stacking interaction potential between adjacent base pairs with a y— dependent stiffness 
constant K. The Morse potential (second term) describes the interaction of the base pairs 
in the two complementary strands. The parameters of the model, the depth D„ and the 
spatial range a„, distinguish between the A-T and G-C base pairs. The integral Eq.(l) can 
be evaluated in the thermodynamic limit using the eigen values and eigen functions of a 
transfer integral (TI) operator [36] 

^ j dyn-i e-^^(^-^"-) </.(y„_i) = e-^^ (3) 



It yields Zy — exp{—(3Neo) where eo is the lowest eigen value of a Schrodinger type 
equation 

1 d'(P{y) 



+ V{y)ct>iy)^ie + so)cl>iy), (4) 



where 



As expected, the constant a does not affect the final expression for Zy or the free energy 
per site / = —{l/N/3) ln{Zy) — eq. However, as will be shown later, it is an important and 
useful parameter of the present formalism which affects the thermal activation energy Sq and 
also it is required to make the TI equation (3) dimensionally consistent. It may be mentioned 
here that the full partition function that contains all the factors including the kinetic terms 
is dimensionless due to the phase space factor h?^ in the denominator. Therefore, the factor 

in (1) can be considered as the part of the phase space distance associated with only 
the configurational part of the partition function whose value can be fixed by comparing 
the results with the experimental measurements. The factor K is assumed of the form 
K — UqII + p{y)] so that K decreases from kQ{l + po) to Uq at large y. It is now possible 
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to separate sq into two parts: a y— dependent positive term which gives an entropic barrier 
when added to the Morse potential and a independent thermal part yth, although both 
parts depend on temperature. Accordingly, Eq.(4) can be rearranged to give 

1 d'<l>{y) 



+ Veff{y)(l){y) = {e + eth) Hy), (6) 



where 

Veffiy) 

and 



D{e-^y-iy-D\+^ln[l + p{y)], (7) 



In the PB model, a y-dependent stiffness with an exponential form p — pQe~^ was 
used to reflect the fact that the DNA is signiflcantly more rigid in the double-stranded 
conformation than when either of the two interacting base pairs is stretched. Although 
there is no theoretical guidehne, a value of 6 = 0.35 A°~^ was used in [11] to simulate the 
melting curves of short heterogeneous DNA sequences ( which accounts for total yield /). 
However, in this study, we have felt the need to use a very small value of h in order to 
explain the spontaneous melting fraction of the fully dissociated molecules and partially 
open intermediate states simultaneously for many oligomers. A small value of h would mean 
a flat entropic barrier over a large distance. Although it is very much essential for p to 
decrease from po to 0, the rate of decrease is not found fast enough unless the stretching 
exceeds a particular limit which may be due to the strong backbone support it receives from 
the neighboring base pairs while in the double stranded conformation. 
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FIG. 1. The effective potential Veff (in units of meV ) as a function of y. The parameters 
are D = 0.05 eV, a = 4.2 T = 300° K, ko = 0.013 eV/A"'^ and po = 2.0. For exponential 

dependence, b is taken as 0.2 A°~^ whereas for Wood-Saxon form, yo = 15 A" and b = 1.42 A°. 

This could be interpreted as some type of co-operative phenomena involving several base 
pairs with a range larger than what is normally expected. To simulate this effect, we use a 
slightly different form of Wood-Saxon type 

Po 



P 



(9) 



1 _|_ Q-{y-yo)/b ' 

characterized by the depth pq, range z/q and diffuseness parameter b. These are the co- 
operativity parameters and, in general, may depend on the sequence composition, length 
of the sequence and also on the buffer conditions. Since it is not known how to fix, these 
are considered as free parameters to be fixed by the experimental observations. Figure 1 
shows a typical plot of the effective potential as a function of ?/ at a temperature of 

T = 300°ir. The potential parameters are D = 0.05 eV, a = 4.2 A°~\ ko = 0.013 eV/A°^ 
and Po = 2.0. For exponential dependence b is taken as 0.2 A°^'^ whereas for Wood-Saxon 
form, uq and b are chosen as 15 A° and 1.42 A° respectively (as an example). The solid 
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line corresponds to the case when the y-dependence is ignored. As can be seen from the 
figure, although both the exponential and the Wood-Saxon forms lead to the same physical 
consequence, namely, a decrease of the backbone stiffness from A;o(l-|-po) to ko, qualitatively, 
both give rise to two different forms of entropic barriers. We have studied both the forms 
and prefer to retain the Wood-Saxon form which is found to be more suitable to explain the 
melting behavior of several oligomers. 

III. DOUBLE STRAND SEPARATION VERSUS INTERMEDIATE STATES 

In the following, we describe the method to characterize the fraction of fully open 
molecules p and partially melted intermediate states a. A dsDNA corformation is ther- 
modynamically stable when the free energy per base pair e < 0. On the otherhand, for 
ssDNA, e > 0. The bound state solution of the Morse potential with the entropic con- 
tribution with inclusion of both harmonic and anharmonic stacking interaction has been 
studied in the past. For long DNA sequence, the denaturation begins with formation of 
bubbles due to local meltings of different domains (local openings of base pairs) where the 
base pair bindings are relatively weak. When the average bubble size attains a critical 
value at temperature close to the melting point (T^), the entropy gain becomes sufficient 
to overcome the binding and the entire sequence separates into two single strands discon- 
tinuously. Therefore, it is only at T >> T^, a real separation is expected. However, for 
shorter sequences as in the case of many synthetic oligomers, the complete strand separation 
may occur even below the melting temperature and the melting process may proceed with 
the co-existence of both fully dissociated molecules (p) and partially melted intermediate 
states (a) over a wide range of temperatures [33,34]. It is noticed that the short oligomers 
(< 10 bp) melt suddenly and the melting curve can be represented by a two state model 
in which the molecules are either completely closed or completely open. As the oligomer 
size increases, zipper model seems to be more appropriate which opens continuously giving 
rise to partially open intermediate states [37] . Many synthetic oligomers which are clamped 
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at both the ends by the relatively GC rich base pairs may nucleate from the middle AT 
region continuosly until a limiting bubble size is reached after which the melting becomes 
discontinuos or even depending on the sequence composition, the melting may have contri- 
butions from both the processes simultaneously. It is also noticed that the oligomers which 
are open at one end behave differently than when both the end points are fixed. Therefore, 
the new experimental measurements with different oligomer designs now impose additional 
constraints on the theoretical models which are required to explain both p and a (hence /) 
simultaneously. Moreover, a quantitative understanding of the melting process of different 
synthetic oligomers will be a stepforward in gaining more insight of the process by which 
the base pairs open up locally in a real DNA sequence even at room temperature. In this 
work, we propose thermally assisted barrier penetration phenomena as a new mechanism for 
complete strand separation whereas the partially opened intermediate states are interpreted 
as the bound state solutions of the same effective potential that generates a barrier. These 
are two distinct processes below the melting temperature although this distinction vanishes 
as the temperature is raised up to a limit above which the bound state no longer exists. 
Below the melting temperature, small amplitude coherent fluctuations which involves sev- 
eral base pairs may drive a dsDNA out of the potential well into the continuum for which 
e > 0. Although thermodynamically favorable to remain separated, these conformations are 
still in the dsDNA state (small amplitude fluctuations around the equilibrium value y ~ 0) 
and are required to overcome an entropic barrier (to the right) in order to remain as stable 
ssDNA conformation. The fraction p, thus, amounts to calculating the thermally assisted 
tunneling probability through an entropic barrier with energy e^h as given by (8) where we 
have set the free energy e = 0. The probability of tunneling can be obtained by solving the 
Schrodinger Eq.(6) using WKB approximation 




(10) 



with the action integral 




(11) 
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Note that the term 1/2P'^ K{y) is the quantum analog of When the thermal energy 

Sth becomes greater than the barrier height Vq, the potential barrier can be approximated 
by an inverted parabola, allowing the use of the analytic expression of Hill and Wheeler 

2ti 



1 + exp{ — {Vo - eth)} 

Oi) 



P 

with 



(12) 



1/2 

(13) 



evaluated at the barrier position y^. 

Now, using (11) for e <Vo and (12) for e > Vq, the ssDNA fraction p as a function of T 
can be estimated. For the intermediate states, we estimate the quantity 



/ \My)?dy, (14) 



N 

where is a normalization constant and 0o(y) is the ground state eigen function of the 
effective potential V^// for temperature T < T^. The integral is bounded by two limits whose 
values are very crucial to reproduce the experimental observation of the peak to valley ratio 
of the partially melted intermediate states [34]. In general, ymin > so that the nucleating 
bubble gains sufficient entropy to grow whereas ymax defines a range upto which the entropic 
barrier exists. It may be mentioned here that we do not include the contribution of (f)o{y) 
beyond ymax although (poiv) is still nonvanishing particularly when T approaches T^. At 
such high temperature, the ground state energy becomes quite close to the continuum and 
the melting may proceed with opening of all the remaining base pairs suddenly. Therefore, 
the fraction of 4>o{y) beyond ymax only contributes to the sudden opening which has been 
already included in Eq.(12) for eth ^ ^o- As will be shown later, the experimental data 
supports a long entropic barrier with range yo ~ Umax- Since we are interested for (f)o{y) 
in the range ymin < ?/ < Umax, where the entropic barrier remains practically constant, the 
ground state solution of the Morse potential with a constant stiffness k = ko{l + po) 

(2(i)'^-V2 _ X 

My) = ^ [r(2d- l)]V2 ^'^^^~^^ "^)ea;p[(d - -)ay]. (15) 

11 



can be used to estimate a (as a good approximation). 

It is well known that a sharp melting transition can not be understood only on the basis 
of the ground state solution of Morse potential with a constant stiffness unless a y-dependent 
stacking interaction is included. In the present formalism, the above desired effect is already 
included in the total yield / which is now the sum of a obtained with a constant stiffness 
as part of the bound state solution and p which is the tunneling probability through a 
dynamically generated entropic barrier that includes stiffness variation. 
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FIG. 2. The fraction of open molecules j? as a function of temperature T (in ^C). The param- 
eters are D = 0.0494 eV , a = 4.2 A°'^ , ko = 0.013 eV/A°^, po = 1.99, yo = 40 b = 4.6 A° 
and a = 2.18 A° for L60S36 sequence and D = 0.0505 eV, a = 4.2 A°-^, ko = 0.013 eV/A"'^, 
Po = 2.01, yo = 15 AO, b = 1.42 and a = 2.18 A° for LA2B18 sequence. The filled circles and 
squares are experimental data points 



Now using the above formalism, we have estimated both the fraction p and a for two 
synthetic oligomers L60B36 and L42B18 which are studied experimentally in ref [34]. Both 
the oligomers are clamped at the ends with identical GC pairs and with an AT rich middle 
region of variable lengths (L60B36: length L= 60 bases, bubble forming region B=36 bases. 
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L42B18: L=42 bases, B=18 bases). Figures 2 and 3 show the estimated fractions of com- 
pletely separated molecules (p) and partially melted intermediate states a as a function of 
temperature T. The parameters used in the calculations are listed in the figure captions. 
The agreement of the theoretical predictions with the experimental measurements is indeed 
excellent. 
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FIG. 3. The quantity a = f — p as a function of temperature T (in "C). The parameters are 
D = 0.0494 eV, a = 4.2 A°-^, ko = 0.013 eV/A"'^, po = 1.99, y^m = 9 and ymax = 40 A^ for 
L60B36 sequence and D = 0.0505 eV, a = 4.2 ko = 0.013 eV/A""^, po = 2.01, y™„ = 9 A° 

and Umax — 

15 A° for L42B18 sequence. The normalization constant is taken as unity. The 
filled circles and squares are experimental data points. 

Although, initially there are too many free parameters, most of them have been fixed 
based on different physical considerations. For example, it is noticed that the tunneling 
probability p is not very sensitive to the Morse potential as long as they are reasonable, but 
depends sensitively on the factor K{y) that decides the nature of the entropic barrier. On 
the otherhand, the fraction a and the peak position have strong dependence on the depth 
and range of the Morse potential. As will be shown later,we fix the Morse parameters that 
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give optimum fits to the a measurements. Accordingly, we take D — 0.0494 eV for L60B36 
and D — 0.0505 eV for L42B18 whereas a is fixed at 4.2 A°~^ for both the sequences. In 
fact, these values are nearly identical to what has been used in [11,32] corresponding to an 
AT base pair. As p is the tunneling probability, there exists a temperature Td at which the 
thermal activated energy eth becomes equal to the barrier height Vq (tunneling probability 
becomes 0.5). Since the contribution of the Morse potential is negligible outside the well, 
it is basically the flat portion of the stacking interaction 1/(2/3) Zn(l + po) that decides the 
height of the entropic barrier. Equating eth with Vq at and using the approximation 
K)~/n(l + po)/(2/3), weget 

T, = ^iMi±M. (16) 

Thus, the above relation puts the constraint on the choice of A; = koil + po) and a. 
For same T^, we may have different combinations of k and a. In practice, k is really not a 
free parameter which can be varied randomly as it also controls the temperature dependent 
effective mass that appears in Eq.(6). Although, we have tried various combinations, it is 
found that the fit is optimum when a is close to 2A°. This is an interesting observation as 
this distance also corresponds to an average stretching beyond which the hydrogen bonds 
start breaking. As we have discussed before, a is a phase space factor associated with the 
configurational partition function. It is probably meaningful to expect a of the order of 
2A° for the configurational changes to occur. Subsequently, we fix a at 2.18^4° which is an 
optimum choice for all the sequences. With a fixed, Eq.(16) now gives a very simple but 
an important relation which can be used to estimate the stiffness parameter k — kQ{l-\- po) 
directly by knowing which can be measured experimentally. For example, from the 
experimental curves, we estimate ~ 68.1°C for L60B36 and T^, ~ 70.36''C for L42B18. 
With the choice of a = 2.18^1°, we get k ~ 0.0389 eV/A"^ and k ~ 0.391 eV/A"^ for the 
above two sequences. Although, knowing the full stiffness k is itself important, yet we do 
not know how to split it between ko and po- However, taking the guidance from ref [11] 
where po was taken ~ 2, we find ka should be ~ 0.013 eV/ A""^ . Subsequently, we fix ko at 
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0.013 eV/A°'^ for all the sequences and consider the co-operative parameter po as a variable 
which, although should be close to 2, can be used as tuning parameter so as to reproduce 
the correct mid point temperature T^. Thus, with a, a, Uq and po fixed, we are now 
left with only two free parameters yo and h which are finally adjusted to reproduce the 
experimental measurements. As listed in the figure captions, the strength po turns out to be 
1.99 for L60B36 whereas it is 2.01 for L42B18. This slight variation is expected as both the 
sequences have different dissociation temperature T^. Similarly, the range and diffuseness 
parameters are found to be 40^4° and 4.6^4° for L60B36 sequence. The corresponding values 
for L42B18 are 1574" and 1.42^4° respectively. Although it is difficult to draw any definite 
conclusion, the length yo seems to depend strongly on the size of the bubble forming AT rich 
region or more precisely on the ratio B/L which are 0.6 and 0.43 respectively for the above 
two oligomers. It may be mentioned here that this observation is more meaningful for these 
two oligomers which have identical endpoints and may not be valid in general. However, 
all sequences studied in this work have a nearly flat entropic barrier over a larger range. 
This effect, probably, can be interpreted as the reminiscent of a long range co-operative 
phenomena. 

The parameters for a are also flxed in a similar way. Since the entropic barrier turns 
out to be nearly constant up to the range yo, c has been estimated using Eq.(14) with the 
ground state wave function given by Eq.(15). It is noticed that the peak position of the yield 
is very sensitive to the melting temperature = 2y/2kD/(kpa) ^. For a given k which is 
fixed from the p yield, \/D/a ratio (the p yield is not sensitive to this ratio), is adjusted to 
have the correct peak position whereas the limit of integrations and i/max are adjusted 
to reproduce the observed peak to valley ratio. As discussed before, we take a — 4.2 A°'^ as 



Note the distinction between and where is the dissociation temperature at which p 
becomes half whereas is the melting temperature above which the bound state vanishes. In 
general, and need not coincide. 
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in [11] and use D = 0.0494 eV for L60B36 and D = 0.0505 eV for L42B18 to reproduce the 
a curves. It is also noticed that the optimum results are obtained when Umax ~ Uo, the 
range of the entropic barrier. Subsequently, we set Umax — Ho and adjust which turns 
out to be 9A° for both the sequences. Therefore, by putting the constraint to explain both p 
and a measurements simultaneously, we have fixed most of the free parameters of the model. 
Although, i/min is expected to be > 0, a value of 9^4° is indeed too large as compared to the 
minimum threshold of average stretching < y >— 2A° above which hydrogen bond between 
the base pairs starts breaking. This large value of ymin probably is very specific to these two 
oligomers whose end points are fixed by a relatively stronger GC rich regions. An unzipping 
transition leading to the continuous melting may begin with a nucleating bubble initiating 
from the softer AT rich middle region. Since this melting is a competition between the energy 
of the bound states and the entropy of the dissociated states, the transition may not proceed 
until the bubble radius attains a minimum value which is about 9A° for these oligomers. A 
stretching distance of ymin ~ would mean a nucleating bubble of diameter d ~ 18^4° 
containing about 4 to 5 open base pairs. Unless the size attains this minimum value, the 
entropy gain is not sufficient enough for the bubble to grow. Note that this minimum size 
is different from the critical bubble size which is required to drive a spontaneous transition. 
The ymin value strongly depends on the size and composition of the sequence. We have 
found that y^in value is much smaller if the sequence is homogeneous. Since the bubble 
forming AT region of L60B36 and L42B18 ohgomers are bounded by relatively stronger GC 
rich regions at the end, the denaturatcd bubble may encounter higher resistance to grow 
due to strong binding. This is to say that the interaction of the bubble with itself and also 
with its surroundings may bring a reduction in entropy. Therefore, for a given entropy gain, 
the minimum bubble size required to grow (by compensating the binding energy) for the 
sequence with fixed end points is larger then the case when the sequence is open at both the 
ends. This could be viewed as analogous to the excluded volume interaction whose effect is 
to reduce the entropy of the loop. Notice that when ymin is large or the range {ymax — Umin) 
is small, the a contribution decreases. Another parameter which may also suppress a yield 
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is the strong binding effect of the base pair. If the binding energy is large, i.e. if the \fD/a 
ratio is large, the temperature at which the bound state vanishes goes up where as the 
dissociation temperature remains unaffected. If » T^i complete dissociation may 
take place much earlier, thus reducing the fraction of continuous melting. Therefore, in 
either case, the most probable route available for denaturation is the sudden melting thus 
making the transition to appear discontinuous. 

It may be mentioned here that for a complete strand separation, it is not necessary to 
raise the temperature until all the base pairs break. Normally, the nucleation is initiated 
with a small denaturating bubble from the softer AT rich region. When the average bubble 
size attains a critical value, the loop entropy becomes sufficient enough to overcome the 
binding so that rest of the base pairs unwind instantly. The critical sized bubble is formed 
either through a sequential unzipping process what is known as a continuous melting or 
through the sudden opening due to coherent fluctuation (followed by tunneling) . In the first 
approach, the critical size is reached when the temperature is close to the melting point 
whereas the second process can happen at any temperature with different probabilities. 
Since the tunneling probability p depends only on the entropic barrier parameters, the 
heterogeneity of the sequence which introduces a site dependent binding Ki, does not affect 
the full dissociation yield. On the otherhand, the fraction of the continuous yield a is 
strongly affected by the heterogeneity of the sequence. By design, the bubble formation for 
L60B36 and L42B18 is initiated in the AT rich middle region. Since the AT rich region for 
L60B36 is long enough, the critical bubble size may still remain confined to the homogeneous 
middle portion only. For L42B18, the homogeneous AT region may not be long enough and 
the bubble may extend upto the GC rich region also. As a consequence, a transition may 
encounter different melting regions with appropriate weights. As can be seen, the a peak for 
L60B36 is more sharper as compared to the case of L42B18 sequence where the GC region 
may start contributing. This curve can be better understood as a weighted sum over two 
melting regions with appropriate weights although p curve can still be understood as a single 
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tunneling phenomenon through an appropriate entropic barrier. 

Figures 4 and 5 show the comparison of /, p and a with the corresponding experimental 
results using the set of parameters as listed in the figure captions 2 and 3. The linear rise 
of the total open base pairs / after the sigmoidal transition is a well-known phenomenon 
attributed to the residual base unstacking in the single strands. 
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FIG. 4. The quantities /, p and o" as a function of temperature T (in "C) for L60-B36 sequence. 
The filled triangles, circles and squares are the corresponding experimental data points. 

However, the theoretical prediction of / does not increase beyond unity as this effect has 
not been included in the present model. 

As shown in figures (2-5), the present model has been quite successful in explaining both 
qualitatively and quantitatively experimental observations. We have also studied two more 
oligomers L42V1 and L48AS for which experimental data points are available [37]. The 
L42V1 sequence is 42 base pairs long with 18 AT base pairs in the middle whereas L48AS 
is 48 base pairs long with a AT rich region of 21 base pairs which is open at one end. 

Figure 6 shows the plot of p as a function of T for the above two oligomers with the 
parameters listed in the figure captions. The parameter is found out to be 18 A° for L42V1 





L60B36 








P 







-f 




.» 








t 






w 

» ' 


■ 






* 


', ■■■■ 


*r , , rr'-TTr. rr.'r. . . . rr-rrvr:-.—:- 



18 



with a diffuseness value of 2.2 A°. This ohgomer {a, B/L ratio of ~ 0.42) is quite identical 
to L42B18 oligomer in design. Accordingly, the parameters of p{y) follows the systematic 
and also the full dissociation probability is quite close to that of L42B18 oligomer. 
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FIG. 5. The quantities /, p and cr as a function of temperature T (in ^C) for 1/42518 sequence. 
The filled triangles, circles and squares are the corresponding experimental data points. 

Interestingly, this oligomer has no contribution to a as it melts suddenly at all temper- 
atures. As mentioned before, the vanishing of a component would mean a vanishing range 
{.Umax — Umin) or a higher effective binding energy which pushes the melting temperature 
Tm to the higher side. Since L42V1 is quite similar to L42B18 in design, it is difficult to 
understand why cr vanishes for the former whereas this component is still significant for 
the later oligomer. Although this needs more careful analyses, a possible reason could be 
different buffer conditions in which the DNA solutions were prepared. It is known that a 
stronger ionic concentration may increase the binding (higher hence higher T^) while the 
entropic parameters may still remain unchanged. On the otherhand, the melting behavior 
of L48AS is quite different from the rest of the oligomers as it opens from one end. It has a 
range which is relatively larger (?/o ~ 60 A° and 6 = 9 A°). The large entropic range explains 
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why p component is small with a relatively larger a contribution. We have not estimated 
the a fraction, but it is expected to be broad with contributions coming from two different 
melting regions. 
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FIG. 6. The fraction of open molecules p as a function of temperature T (in ^C). The param- 
eters are D = 0.0505 eV , a = 4.2 kQ = 0.013 eV/A°'^, po = 2.01, yo = 18 A^, b = 2.2 and 
a = 2.18 A" for L42V1 sequence and D = 0.0505 eV, a = 4.2 A"-'^, ko = 0.013 eV/A"'^, po = 2.07, 
Uo = 60 A^, b = 9 A^ and a = 2.18 A° for LA8AS sequence. The filled circles and squares are 
experimental data points 

The temperature driven denaturation has also been studied for synthetic ploy (A)- ploy 
(T) oligomers of various lengths in different buffer conditions [37]. These oligomers can 
be regarded as homogeneous duplexes held together by identical base pair interaction at 
the different sites. Although, the measurements do not separate the fractions into fully 
dissociated and partially open components, the total fraction / is itself interesting as it has 
a large not so well understood premelting region. The peak and width of the derivative with 
respect to temperature also show strong dependence on the length of the oligomer as well as 
on the buffer conditions. Due to finite size effects, the co-existence of both the components 
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is expected for these oligomers also. Since data is not available, we will not try to estimate 
p and a exactly so as to reproduce the total yields, rather we will discuss a scenario based 
on which the experimental results can be understood at least qualitatively. Again as before, 
we consider an arbitrary oligomer which may have similar potential parameters like L60B36 
oligomer except for the fact that this oligomer is now homogeneous. A homogeneous sequence 
may have a lower ymin value than when the sequence composition is more heterogeneous. 
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FIG. 7. The quantities /, p and cr as a function of temperature T (in ^C) for an arbitrary 
sequence. The parameters are D = 0.0494 eV, a = 4.2 ko = 0.013 eV/A°'^, po = 1.99, 

2/0 = 40 A^, b = 4.6 ^0, Q = 2.18 A°, y^in = 0, Vmax = Vo and N = 0.5. 

As a model study, we choose i/min = 0, i.e. the minimum expected value. As the i/min 
decreases, the behavior of a changes completely, becoming more uniform as ymin decreases 
from a higher threshold value towards 0. Figure 7 shows a typical plots of a, p and / 
as a function of T where we have set ymin to its lowest value. Although we have used a 
normalization factor N = 0.5 , it can be adjusted so as to reproduce the appropriate yield 
in the premelting region. The total yield is quite similar (at least qualitatively) to what 
one observes experimentally [37]. Interestingly, the so called premelting region is probably 
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due to a continuous transition which has nearly uniform probabihty up to temperature 
above which the transition becomes discontinuos. Since below T^, the continuous process 
dominates, a zipper type model may be found more appropriate to describe this type of 
melting curve, particularly when the sequence is homogeneous and also relatively large. It is 
also observed that the length of the entropic range yo depends on the length of the sequence, 
i.e. for the same composition, yo increases with the length of the sequence probably until 
an upper limit is reached. So the smaller sequence will have short ranged entropic barrier 
with increased p contribution. When the size is very short (< 10 bp), it may completely 
dissociate without any intermediate state. This explain why a two step model is found more 
appropriate to describe the melting behavior of very short oligomers. 

Since the premelting region is fiat, the non vanishing part of the derivative of the total 
yield is basically due to the fraction p. Figure 8 shows the plot of dp/dT (estimated numer- 
ically) as a function of T for different co-operativity parameters. The curves on then left 
correspond to the case when po — 1-8 whereas the curves on the right are for po — 1-99. On 
both the sides, the solid curve corresponds to the parameters as listed in the figure captions 
while the dotted and dashed curves are obtained when either yo or ^ is changed. The point 
which needs to be stressed here is that the peak positions (hence the temperature T^) are 
quite sensitive to the strength po which may be small for smaller oligomers. Therefore, the 
shift in peak positions of the oligomers to the higher values when the length increases can be 
understood on the basis of the increase of the po parameter. On the otherhand, the width 
and the height are very sensitive to t/o ^-^d h parameters which also depend on the length 
of the oligomers as well as on the buffer conditions. Therefore, with suitable po, yo ^-nd h 
parameters, it is possible to understand all aspects of the melting behavior of the synthetic 
oligomers studied in [37]. 
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FIG. 8. The derivative of p as a function of temperature T (in ^C) for some arbitrary sequence. 
The parameters are D = 0.0494 eV , a = 4.2 Uq = 0.013 eV/A""^, a = 2.18 yo = 40 A° 

and h = 4.6 The curves on the left correspond to a po value of 1.8 whereas the curves on the 
right are for po = 1.99. In both sets of curves, the solid lines correspond to the parameters as listed 
above. The dotted lines are obtained when only hq is changed from 40 A° to 32 A° and the dashed 
lines are obtained when h is changed from 4.6 A" to 3.5 A°. 



Based on the above studies, we have several conclusions to make. Recall again, a large 
DNA sequence melts with the initiation of the nucleating bubbles at different regions which 
are relatively soft. A complete strand separation occurs when the average bubble size attains 
a critical value. Depending on the nature of the entropic barrier, the domain melting occurs 
either through continuous opening of base pairs or due to coherent fluctuation followed by 
tunneling. However, at this level, whatever may be the process of domain melting, the 
complete dissociation of the full sequence will not be realized unless the average bubble size 
attains a limiting value so that the entropic gain becomes sufficient to overcome the bound 
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state energy. Once the melting point is reached, a complete strand separation occurs what 
appears as a discontinuous phase transition. Therefore, the complete DNA denaturation in 
the thermodynamic limit can be interpreted as a first order phase transition characterized by 
a sharp melting temperature at which the bound state vanishes. However, due to finite size 
effect as in the case of oligomers, both type of processes, one due to continuous opening and 
the other due to the sudden melting, may co-exist and the concept of a phase transition is 
not strictly applicable. Therefore, we invoke a barrier penetration phenomena along with the 
bound state solution to explain these two fractions separately. Although, it is not possible 
to draw any definite conclusions unless more experimental data sets of various synthetic 
oiligomers are available with different design, we can only make a few qualitative remarks 
as follows. The available measurements seem to support a long ranged entropic barrier 
with parameters which depend on the length and composition of the sequence. For the 
same length, the range of the sequence made up of GC base pairs is higher than when the 
sequence is homogeneous with AT base pairs. Consequently, the GC rich sequence will have 
less p and more a contributions as compared to the AT rich sequence. If the oligomer is 
short enough (< 10 bp) and homogeneous, a sudden melting may occur due to an entropic 
barrier which is very short ranged (small value of yo). These oligomers will also not have 
any continuous fraction due to vanishing {ymax — Vmin) range and the full melting process 
can be explained by the so called two state model. On the otherhand, with increasing size 
of the oligomers, the entropic barrier range increases (due to stronger collectivity) and the 
full dissociation probability p decreases. If ymin is close to zero, continuous opening takes 
place with equal probability at all temperatures until the melting point is reached. As 
seen in [37], these type of processes are characterized by a slowly rising pre-melting region 
followed by a sudden transition and are more suitable to be described by the zipper type 
model as contrast to the two state phenomena. The melting behavior of the oligomers 
studied in [34] can also be understood as the sum of two components. These two oligomers 
L60B36 and L42B18 with the B/ L ratios of 0.6 and 0.43 are clamped at both the ends with 
identical GC rich regions. Since the relative AT length is large enough, the bubble forming 
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region is homogeneous and the corresponding yo values of 40A° and 15A° are probably 
following the B/L (relative length) systematic. The L60B36 has a longer entropic barrier 
than L42B18 and the corresponding fraction p is also less for L60B36. On the otherhand, 
the entropic range is smaller for L42B18. Accordingly, the a contribution is also less as 
compared to L60B36. However, the key aspect where these two oligomers differ from not 
being homogeneous is the higher ymin value of about 9^4° which makes the peak to valley 
ratio in the a yield very distinct. The oligomer L42V1 which is nearly identical to L42B18 
has also been studied in [33]. While this oligomer has a comparable p yield, it does not 
have any a component. This aspect is interesting and needs further investigations. On the 
otherhand, the ohgomer L48AS has a, B/L ratio of 0.25 with more GC components and 
also with a large entropic range. Consequently, the p component for this type of oligomer 
is negligible and the transition mostly proceeds through continuous opening. For the same 
reason, as remarked in [34], even shorter oligomer (down to 19 bp) which unzips from one 
end may open continuously. 

Finally, we would like to add here that the shape of the intermediate fraction yield a 
depends strongly on the sequence heterogeneity since becomes site dependent. However, 
the full dissociation probability p is still can be obtained from the barrier penetration for- 
malism as it is sensitive only to the entropic barrier. The temperature at which the full 
dissociation probabihty becomes half is directly related to the DNA stiffness parameter. For 
the sequence for which — T^, the mid point of the full absorption curve / can be used 
to estimate the stiffness parameter. In general, need not be same as T^. Since relates 
directly to the stiffness parameter and also its derivative is more sensitive to the parameters 
of the entropic barrier, it will be more meaningful to study these two components p and a 
separately. 
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APPENDIX: A 



A linear chain of DNA molecules can be described by the Hamiltonian given by 

^ = E l^i'^l + ^n) + 1^ [(«n - Un-l)^ + {v^ - V^-lf] + V{Un - V^), (M) 



where Un and v„ correspond to the displacements of the bases at site n along the direction 
of the hydrogen bond, and V is the potential between the base pairs which can be of the 
Morse type 



V{Un -Vn)^D [exp {-A{Un - Vn)} - if ■ 



(A2) 



In terms of the variables Xn = (wn + Vn)/V^ and yn = {un — Vn)/V^, the above Hamil- 
tonian can be written as 



where Pn = mXn and g„ = mi/n are the canonical momenta, and 
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(A3) 



(A4) 



(A5) 



The statistical mechanics of the model is described by the partition function 
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1 

^^¥^J-oo ^^=^ d^ndyndpndqn e-^^(^'"'-"'«"'^"), (A6) 

where P = l/{kpT) with kf^ being Boltzmann's constant, and h is Plank's constant. The 
integrals for the variables p„, g^, and Xn are of the Gaussian type, and give, respectively 

Zp = Z,= {27rmkf,Tf/' , = {2T:kpT/Kf'^ . (A7) 

The remaining configurational partition function Zy involves nonlinear Morse potential, and 
given by 

Zy^^ Ug, dyn e-^^(y-y--^\ (A8) 
a" J-oo 



The factor a (dimension of length) has been in introduced to make Zy dimensionless. 
The above integral can be evaluated by using the transfer integral (TI) method defined by 

- / dvn-i e-^^(^-^"-) = e-^^ <P{yn). (A9) 

a J 



To compute the integral, we consider a finite chain of N base pairs and chose periodic 
boundary conditions by demanding that yjsf — yo- This amounts to adding a fictitious 
base pair with index which has the same dynamics as base pair N, with a simultaneous 
introduction of a factor S{yN — yo)- Consequently, Eq.(A8) can be written as 



— ^ J ^n=0 dyn 



n 



-PH{yn,y„-i) 



HyN-yo)- 



(AlO) 



The 5 function of Zy can be expanded on the base of the eigenfunctions of the TI 
operator as S{yN - yo) = Ei 4>i{yN) (piiyo)- Performing the integrals over yo, yi, y2,....,yN-i 
successively yields 

Zy ^ i: e-'^^- / dy^ \UyN)\' = E^"'"''^'- (All) 

i i 
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In the thermodynamic hmit {N — > oo), the result is dominated by the term having the 
smallest value of ej which we denote by eo- The free energy per site is therefore 

/ = -^ln{Z,) = 60. (A12) 

More interesting for the study of DNA denaturation is the mean stretching < y > of the 
hydrogen bonds. It is given by 

<y> = Yj^n=iy e-^^^^y-'^-^Uy^. (A13) 
The integral can again be calculated with TI method and yields 

^ ^ _ Eili <Uy)\y\Hy)> e'^^-^ 

^ Ef=i <My)\Uy)> e-^^- ■ ^""'^^ 

Again, in the limit of large N and considering only the lowest eigen value eo, the average 
stretching is given by 

<y> = <My)\y\My)> ^ j \(t>o{y)?dy, (Ais) 

for the normalized eigenfunction (poiy)- Note that with the transformation of the co- 
ordinates, < y > gives the mean stretching of the hydrogen bond between the two 
opposite base pairs. 

The problem is therefore to find the eigenfunctions and eigen values of the TI operator. 
The calculation can be made analytically in the limit of strong coupling between sites. For 
large K, when y„ differs significantly from yn-i, the term K{yn — yn-if grows rapidly, so 
that the TI is dominated by the value of which are close to y^. By defining yn-i — yn+z 
and expanding (j){yn + z) in powers of z, Eq.(A9) can be written as 



1 r°° ^Kz^ 
e 2 



a J-c 



■ d<l>{y) 1 ,d'<j,{y) 

(j){y) + z— — + -z 



^-Ple-ViV2y)] (A16) 



dy 2 dy^ 

where the index n has been dropped. Since the odd term of the Gaussian integral in z 

vanishes, the above integral reduces to 

1 



/ 27r 



2pK dy^ 
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The factor ^^^2 can be absorbed in the eigen value by redefining 



In many practical situations, the relevant magnitude oi {e — V) ^ D , the depth of the 
Morse potential. For positive y (as in the case for stretching), the Morse potential is bounded 
by D. If /3D < 1, the exponential can be expanded to get Schrodinger type equation 

1 920(1/) 



where 



+ ny)0(y) = e0(y), (A19) 



V{y)^ D\e-"-y -if , a^V2A, (A20) 



and e as defined by Eq.(Al8). 



The above equation has a discrete spectrum when d = \/2KD / {kpTa) > 1/2 and the 
eigenvalues and the normalized eigen functions of the ground state are then 

(A21) 



2/3 \ 2n J P \2KJ SP'^K' 
(2(i)'^-V2 



My) = V^.^^/-^exp{-de-'^y)exp[-{d - -)ay]. (A22) 
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